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Abstract 

We present a clear A^-body realization of the growth of a Bahcall-Wolf f °c E^''^ (p °^ r^^'^) density cusp 
around a massive object ("black hole") at the center of a stellar system. Our A^-body algorithm incorporates 
a novel implementation of Mikkola-Aarseth chain regularization to handle close interactions between star and 
black hole particles. Forces outside the chain were integrated on a GRAPE-6A/8 special-purpose computer 
with particle numbers up to A^ = 0.25 x 10^. We compare our A^-body results with predictions of the isotropic 
Fokker- Planck equation and verify that the time dependence of the density (both configuration and phase-space) 
predicted by the Fokker-Planck equation is well reproduced by the A^-body algorithm, for various choices of 
A^ and of the black hole mass. Our results demonstrate the feasibility of direct-force integration techniques for 
simulating the evolution of galactic nuclei on relaxation time scales. 
Subject headings: stellar dynamics, galaxies: nuclei, black holes 



1. introduction 

The distribution of stars around a massive black hole is a 
classic problem in galactic dynamics. Many such distribu- 
tions are possible, depending on the initial state, the mode 
of formation of the black hole, and the time since its for- 
mation. However, a number of plausible scenarios predict 
a steeply-rising stellar density within the black hole's sphere 
of influence, r <, rh = GMbhI'^'^, with Mbh the black hole 
mass and a^ the ID stellar velocity dispersion outside of r/,. 
Such models receive support from the observed run of stel- 
lar density with radius near the centers of the nearest galaxies 
known to contain massive black holes: the Milky Way, M31 
and M32. The nucleus of each galaxy has a density cusp with 
p ~ r^^, Y~ 1-5 within the black hole's sphere of influence 
JLauer et al.ll998HGe nzel et a lJ2003l) . This is consistent with 
the slope predicted by so-called adiabatic growth models, in 
which a black hole with small initial mass is embedded in 
a star cluster and its mass increased to some final value, on 
a time scale long compared with orbital periods. The final 
density in the adiabatic growth model follows a power law at 
r<,rh, with index y/ that depends on the initial stellar distribu- 
tion function. If the initial model is a non-singular isothermal 
sphere, then y/ = 3/2 (Peebles 1972a; Young 1980). However 
a non-singular isothermal sphere seems a rather ad hoc guess 
for the initial state, and initial mod els without fla t cores give 
steeper final slopes, 1 .5 ;$ y/ ;$ 2.5 llMerriti2004 . 

Uncertainties in the initial state are less consequential if the 
stellar cluster is old compared with the time Tr for gravita- 
tional scattering to redistribute energy between stars. In this 
case, one expects the collisional transport of mass and energy 
to set up a steady-state distribution whose functional form is 
indendent of the initial phase-space density. Peebles (1972b) 
first addressed this problem and derived a power-law index 
y = 9/4 for the stellar density within r/,. Peebles obtained his 
solution by setting the flux of stars in energy space to a con- 



stant non-zero value. Shapiro & Lightman (1976) and Bah- 
call & Wolf (1976) criticized the Peebles derivation on the 
grounds that the implied flux is unphysically large, in fact di- 
vergent if the solution is extended all the way to the black 
hole. A physically more reasonable solution would have a 
nearly zero flux of stars into the black hole. Bahcall & Wolf 
(1976) repeated Peebles's derivation, setting the phase-space 
density f{E) to zero at the black hole's tidal disruption radius. 
They found that / evolves, in a time only slightly longer than 
T,-, to a steady state in which the flux is close to zero at all en- 
ergies. The zero-flux solution has / ^ E^''^ within the black 
hole's sphere of influence, with E >Q the binding energy; the 
corresponding stellar density is p 0= r^^'^. 

The Bahcall-Wolf solution has been verified in a num- 
ber of subsequent studies, almost all of which wer e 
based on the Fokker-Planck formalism ("Cohn & Kulsru 
fT978t iMarchant & Shapirol (198(1 iFte itag & Benz 200" 
lAmaro-Seoane. Freitag & Spurzerrl2004h . The / 0= £''^ and 
p oc r^il'^ character of the solution has been found to be ro- 
bust, at least at radii where the capture or destruction of stars 
occurs in a time long compared with orbital periods. But ver- 
ifying the Bahcall-Wolf solution in an A^-body integration is 
also clearly desirable, since an A^-body simulation is free of 
many of the simplifying assumptions that limit the applica- 
bility of the Fokker-Planck equation, including the restriction 
to small-angle scattering, and the neglect of spatial inhomo- 
geneities in the derivation of the diffusion coefficients. But the 
A^-body approach is challenging: large particle numbers are 
required to resolve the cusp, and accurate integration schemes 
are needed to follow the motion of stars near the black hole. 

In this paper we combine a sophisticated A^-body code with 
the special-purpose GRAPE hardware and show that the for- 
mation of a Bahcall-Wolf cusp can be convincingly repro- 
duced without any of the approximations that go into the 
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Fokker-Planck formalism. Our results provide a clear demon- 
stration of the applicability of direct-force A^-body techniques 
to the collisional evolution of dense star clusters, and high- 
light the usefulness of direct A^-body techniques for under- 
standing the dynamical evolution of galactic nuclei containing 
supermassive black holes. 

2. METHOD 
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algorithm is an adaptation of NBODYl 
to the GRAPE-6 special-purpose hardware. 
Forces for the majority of the particles were computed via 
a direct-summation scheme without regularization. Particle 
positions were advanced using a fourth-order Hermite inte- 
gration scheme with individual, adaptive, block time steps 
jAarsethll2003h . Time steps of stars near the central massive 
particle ("black hole") would become prohibitively small in 
such a scheme. Therefore, the motion of stars near the black 
hole was treated via the chain regularization algorithm of 
Mikkola and Aarseth (Mikkola & Aarsethl[T9 9Q. 1993). This 
scheme selects an ordering for the n particles in the chain 
and then regularizes the n — 1 pairwise interactions accord- 
ing to the standard Kustaanheimo-Stiefel (KS) transformation 
jKustaanheimo & Stiefel 1965), such that motion in Cartesian 
space with a singular potential is transformed into harmonic 
motion in a 4D vector space. 

Our chain method differs from the standard one (which is 
incorporated in NBODY5 and NBODY6) in two ways. First, 
the massive body was always positioned at the bottom of the 
chain; and second, the chain was ordered according to force, 
not distance. The chain subsystem typically consisted of the 
black hole plus a small number of stars, n — 1 w 10 — 20. The 
chain's center of mass was a pseudo-particle as seen by the 
A^-body code, and was advanced by the Hermite scheme in 
the same way as an ordinary particle. However the struc- 
ture of the chain (i.e. the forces from individual chain mem- 
bers) was resolved as necessary for the accurate integration 
of nearby particles. Furthermore the equations of motion of 
particles within the chain included forces exerted by a list 
of external perturbers. Whether a given star was listed as a 
perturber was determined by a tidal criterion, with typically 
^10^ particles satisfying the criterion. Particles were de- 
noted perturbers if their distance r to the chain center of mass 

satisfied r < i^eriti = ('«/'Wch)'''^Ynim -'^ch. where m^i and/^ch 
represent the mass and spatial size respectively of the chain, 
m is the mass of one star, and y,,,,-,, was chosen to be 10^^; 
thus /?criti ~ 10^(m/m(;h)''^/?ch- The chain was resolved for 
stars within a distance Raia = ^pRch from the chain's center 
of mass, where A,p = 100. In the standard chain implemen- 
tations in NBODY5 and NBODY6 (with particles of nearly 
equal mass), /?ciiti ~ ^ciit2, while in our application these two 
radii are significantly different. 

The Mikkola- Aarseth chain algorithm becomes ill-behaved 
in the case of extreme mass ratios between the particles, 
Mrh Im ^ 10^ (S. Mikkola, private communication). Aarseth 
(|2003) advocates a time-transformed leap frog (TTL) to treat 
such extreme mass ratios. The largest mass ratio in our simu- 
lations is '^ 10^ for which the standard chain algorithm works 
perfectly well. Using a minimum tolerance of 10^'^ for the 
chain's Bulirsch-Stoer integrator allowed us to reach typical 
relative accuracies of 10"^. 

3. INITIAL CONDITIONS 

A crucial element of our method is the use of initial con- 
ditions that represent a precise steady state of the collision- 



less Boltzmann equation. Embedding or growing a massive 
particle in a pre-existing stellar system can easily result in 
the formation of a density cusp with slope similar to that 
predicted by Bahcall and Wolf (1976), but for reasons hav- 
ing nothing to do with collisional relaxation. For instance, 
as discussed above, adiabatic growth of a black hole pro- 
duces a density profile within the hole's sphere of influence 
of p ~ r^'^'p, 1.5 ^ Yvp ^ 2.5; hence a cusp that forms colli- 

sionlessly can easily mimic a p <=<: r^^/"* collisional cusp. To 
avoid any possibility of non-collisional cusp formation in our 
simulations, which would seriously compromise the interpre- 
tation, we generated initial coordinates and velocities from 
the steady-state phase space density f{E) that reproduces the 
Dehnen (1993) density law in the gravitational potential in- 
cluding both the stars and the black hole. The Dehnen model 
density satisfies p(r) tx r^y at small radii, and the f{E) that 
reproduces Dehnen's p(r) in the presence of a central point 
mass is non-negative for all y > 0.5; hence y — 0.5 is the flat- 
test central profile that can be adopted if the initial conditions 
are to represent a precise steady state. As shown in Table 1, 
most of our runs used initial conditions with this minimum 
value of y. Henceforth we adopt units such that the gravita- 
tional constant G, the total stellar mass M, and the Dehnen 
scale length a are equal to one. 

Table 1 also gives the other important parameters of the A^- 
body integrations. The influence radius of the black hole r/, 
was defined as the radius at which the enclosed stellar mass at 
f = was equal to twice the black hole mass. (This is equiv- 
alent to the more standard definition r/, = GMbhI'^^ when 
p o<: r^^.) The relaxation time Tr was computed at f = 0, r = r/, 
from the standard expression (eq. 2-62 of Spitzer 1987), set- 
ting InA = ln(r/,o^/2Gm). This definition of A is equivalent 
to equating b,nax, the maximum impact parameter for encoun- 
ters in Chandrasekhar's theory, with r/,. This choice is mo- 
tivated by the expectation that b„^ax for stars near the center 
of strongly inhomogenous stellar system should be of order 
the radius at which the density falloff begins to abate, e.g. 
the core radius if there is a core (M aoz 1993: Merritt 2001:i) . 
In our simulations, this radius is of order the Dehnen scale 
length fl « a few x r/, at f = 0, decreasing to a fraction of 
r/, after formation of the collisional cusp; hence a choice of 
bmax ~ ri, seems appropriate. We stress that Chandrasekhar's 
InA is a poorly-defined quantity in strongly inhomogeneous 
and evolving systems, and our choice is at best approximate. 
Nevertheless, we will see that the time scaling determined 
by this choice of InA results in very good correspondence 
between the evolution rates seen in the A^-body and Fokker- 
Planck models. 

4. FOKKER-PLANCK MODELS 

We compared the A^-body evolution with the predictions of 
the time-dependent, orbit-averaged, isotropic Fokker-Planck 
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equation: 



3/ 



4n'p{E)^ = - 



dFF. 



3/ 



FEiE,t) = -DEE^-DEf, 



(la) 
(lb) 



DEEiE) = 64iCG''m-'\nA x 
q{E) [ dE'f{E')+ rdE'q{E')f{E') 



De{E) = -64%^G^n?\nA / dE' p{E')f{E') 



(Ic) 



Spitzer 1987). Here p{E) = 
dq/dE is the phase space 



(e.g. Cohn 1980; 

4\/2 j;""" '^' r2 ^E-^{r)dr 
volume accessible per unit of energy, the upper integration 
limit is determined via — 'I>(rmax) = E, and E — — v^/2 — 
<I>(r) > is the binding energy per unit mass (hereinafter the 
"energy.") Since the gravitational potential changes as the 
stellar distribution evolves, equation should also contain 
a term describing the adjustment in / due to changes in E 
as the gravitational potential evolves (e.g. Spitzer 1987, eq. 
2-86). However changes in p(r) only take place well within 
r/, in our simulations, and the potential is dominated by the 
(fixed) mass of the black hole in this region. Hence, assuming 
a fixed potential in the Fokker-Planck equation is an excellent 
approximation for our purposes. 

When scaling the Fokker-Planck results to the A^-body re- 
sults, the only free parameter is In A. We used the values given 
in Table 1 . 

5. RESULTS 

As the cusp develops, the density of stars atr<, r/, increases. 
This is illustrated in Figure 1, which shows the mass in stars 
within a radial distance O.lr/, from the black hole as a func- 
tion of time for each of the runs. (Distances were defined 
with respect to the instantaneous position of the black hole 
particle; the latter wanders like a Brownian particle, but the 
density peak tends to remain centered on the black hole as it 
moves.) For comparison, we also show in Figure 1 the same 
quantity as computed from the Fokker-Planck equation. The 
time scaling of the Fokker-Planck equation was set using the 
value of In A given in Table 1 - no adjustments were made to 
optimize the fit. (We note that integrations like these could in 
principled be used to evaluate InA.) While there are hints of 
systematic differences in some of the runs, overall the corre- 
spondence is very good: clearly, the A^-body evolution is quite 
close to what is predicted from the Fokker-Planck equation. 

The A^-body runs in Figure 1 exhibit a range in A^ from 
0.25 X 10^ (Run 1) to 0.1 x 10*^ (Runs 2-4), showing that the 
correspondence between A^-body and Fokker-Planck results 
remains good over at least a modest range in particle num- 
ber. Figure 1 also suggests that each of the integrations has 
reached an approximate steady state with regard to changes in 
the density at the final time step. 

Having demonstrated the reliability of the time scaling of 
the Fokker-Planck equation, we can make a more detailed 
comparison with the A^-body models. Figure 2 shows the evo- 
lution of the stellar density p (r) in Run 1 compared with the 
Fokker-Planck prediction. The A^-body density was computed 
from snapshots of the particle positions (no time averaging) 
using MAPEL, a maximum penalized likelihood algorithm 
JMen-iU 1994: Memtt & Tremblav 1994) . The radial coor- 
dinate was defined again as distance from the black hole. The 
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Fig. 1 . — Evolution of the mass in stars within a distance O.lr/, from the 
black hole, where r;, is the influence radius of the black hole measured at time 
zero. Noisy curves are from the W-body runs; smooth curves are solutions to 
the Fokker-Planck equation. 




Fig. 2. — Evolution of the mass density profile around the black hole. Left 
panel: W-body; p(r) was estimated from the particle positions at times t = 
100,200,300,400,500. 1500 via a maximum penalized likelihood algorithm. 
Right panel: Densities predicted from the Fokker-Planck equation at the same 
times; scaling of the time unit used the value of InA given in Table 1 . Lower 
dashed curves show the density at f = 0; upper dashed curves show p « r"' ", 
the asymptotic solution to the Fokker-Planck equation. 



correspondence between A^-body and Fokker-Planck results 
is again quite good; the only systematic difference appears at 
very small radii (r < O.Olr/,) where the particle numbers are 
too small for reliable estimates of p. The final cusp is well 
represented by p o= ^^'-^^ ^t r ;=5 O.lr/,, both in Run 1 and in 
the other integrations (not illustrated here). 
The large number of particles in our A^-body experiments 
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Fig. 3. — Evolution of the phase-space density of stars around the black 
hole in Run 5. Left panel: W-body; f{E) was estimated from the particle 
energies at times I = 1 80, 300, 600 as described in the text. Right panel: Den- 
sities predicted from the Fokker-Planck equation at the same times; scahng 
of the time unit used the value of In A given in Table 1 . Lower dashed curves 
show f{E) at ( = 0; upper dashed curves show / « £ ' , the asymptotic so- 
lution to the Fokker-Planck equation. 



permits us to go one step deeper in comparing A^-body with 
Fokker-Planck results. It is possible to extract estimates of 
/(£■) from the A^-body data sets. We did this as follows. Snap- 
shots of the particle positions and velocities were stored at 
each A^-body time unit; such data sets are essentially uncor- 
rected at radii near r/,. Roughly 70 of these snapshots were 
then combined into a single data file, giving an effective A^ 
of ~ 10^. From the combined data set we computed an esti- 
mate of the gravitational potential using standard expressions, 
then computed the phase-space volume element p{E) defined 
above. The particle energies were also computed, and a his- 
togram constructed of N{E). Finally, we used the relation 
f{E) — N{E)/4n'^p{E) to compute an estimate of the phase- 
space density, assuming an isotropic distribution of particle 
velocities. 

Figure 3 shows the results for Run 5 (A^ = 1.5 x 10^, y = 
0.5). While f{E) is clearly harder to estimate than p(r) - 
it is effectively a 3/2 derivative of p and hence noisy - we 
can see from Figure 3 that the / extracted from the A^-body 
runs evolves in a very similar way to the / computed via 
the Fokker-Planck equation, and that its steady-state form is 
consistent (modulo the noise) with the Bahcall-Wolf solution 
/ oc £'1/4 at late times. 



6. DISCUSSION 

Our results contribute to the growing body of literature 
demonstrating the feasibility of direct-force A^-body tech- 
niques for simulating the evolution of dense stellar systems 
on collisional time scales. Two landmark studies in this 
vein were the recovery of gravothermal oscillations using the 
GRAPE-4 special purpose computer by Makino (1996), and 
the computation of the parameters of core collapse by Baum- 
gardt et al. (2003) using the GRAPE-6. Our study extends 
this work to the yet more difficult case of star clusters around 
massive singularities, in which accurate and efficient compu- 
tation of the forces is particularly critical. We have shown 
that particle numbers accessible to a single special-purpose 
computer ( < 10^), combined with a regularization scheme for 
close encounters to the black hole, are sufficient to test predic- 
tions derived from a macroscopic formulation of the evolution 
equations, and that the two approaches indeed make consis- 
tent predictions about the formation of a cusp on relaxation 
time scales. 

This success points the way to even more challeng- 
ing A^-body studies of galactic nuclei. For instance, the 
"loss cone" problem - the rate at which a central black 
hole captures or destroys stars - has so far been the al- 
most exclusive domain of Fokker-Planck treatments. Veri- 
fying the Fokker-Planck predictions via direct A^-body sim- 
ulations is feasible but will require particle numbers large 
enough that loss-cone refilling takes place on a time scale 
long compared with orbital periods. This implies A^ > 
10^ (tMilosavhevic & Merritt 2003), currently achievable via 
distributed hardware (Dorband, Hemsendorf & Merritt 2003). 
Related problems requiring comparably large A^ are the long- 
term evolution of binary supermassive black holes, and the 
Brownian motion of black holes in dense nuclei. 
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